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SECTION  1 
INTRODUCTION 

This  report  surveys  an  emerging  technology  spawned  by  increasing  demands  for  multi¬ 
disciplinary  analysis  capabilities  in  computational  mechanics.  Although  ad  hoc  efforts  in 
mechanics  date  back  many  years  and  systematic  efforts  have  been  made  in  other  disciplines, 
development  of  the  technology  in  mechanics,  herein  termed  “partitioned  analysis  of  coupled 
mechanical  systems”,  only  began  in  about  1975  [l]. 

1.1  COUPLED  SYSTEMS  DEFINED. 

What  do  we  mean  by  coupled  systems?  Within  the  context  of  this  report,  coupled  systems 
are  continuous  mechanical  systems  described  by  partial  differential  field  equations  that 
are  coupled  at  their  boundaries,  which  coupling  is  specified  in  terms  of  nonhomogeneous 
boundary  conditions.  Typical  generic  examples  are: 

•  A  structure  submerged  in  a  solid  or  fluid  medium 

•  Connected  structures 

•  A  porous  solid  permeated  by  a  fluid 

Furthermore,  we  focus  here  on  the  transient  dynamic  interaction  of  coupled  systems,  so 
that  the  governing  equations  for  the  problems  considered  are  sets  of  partial  differential 
equations  in  space  and  time  that  are  coupled  through  partial  differential  boundary  equa¬ 
tions.  Finally,  we  restrict  our  attention  to  situations  in  which  the  field  equations  have 
all  been  semi-discretized,  i.e.,  discrete-element  (finite-difference,  finite-element,  boundary- 
element)  techniques  have  been  used  to  treat  spatial  dependence,  yielding  sets  of  ordinary 
differential  equations  in  time  coupled  by  ordinary  differential  boundary  equations.  We 
refer  to  these  equations  as  the  “discrete-  element  equations  of  motion”.  * 

1.2  SOLUTION  METHODS  FOR  COUPLED  SYSTEMS. 

Four  principal  methods  are  readily  identified  for  solving  the  discrete-element  equations  of 
motion  for  coupled  systems;  they  are: 

•  Explicit  numerical  time  integration  of  the  global  equations 

•  Implicit  numerical  time  integration  of  the  global  equations 

•  Explicit/implicit  numerical  time  integration  of  modal  equations 

•  Partitioned  numerical  time  integration  of  the  global  equations 


*  Limiting  our  attention  to  mechanical  systems  derives  from  the  computational  mechanics 
environment  from  which  partitioned  analysis  has  emerged  and  the  audience  to  which  this 
report  is  addressed;  the  technology  is  readily  extended  to  treat  other  coupled  systems,  how¬ 
ever.  The  focus  on  transient  problems  is  less  restrictive  than  it  might  appear,  in  that  both 
static  and  eigenvalue  problems  may  be  cast  as  evolutionary  ones  (see,  for  example,  [7]).  The 
restriction  to  semi-discretized  problems  is  also  modest,  in  view  of  the  power  and  breadth  of 
discrete-element  analysis  (see,  e.g.,  [22]). 
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Global  Explicit  Integration'.  This  first  method  is  the  simplest,  both  conceptually  and  with 
regard  to  computational  implementation.  We  merely  apply  an  explicit  time-  integration  al¬ 
gorithm  (such  as  the  central-difference  algorithm)  to  the  global  discrete-element  equations 
of  motion  (hereafter  called  the  “global  equations”).  Because  the  response  at  a  given  time 
step  is  extrapolated  from  the  responses  at  previous  steps,  we  can  usually  avoid  the  simul¬ 
taneous  solution  of  many  algebraic  equations,  which  yields  high  computational  efficiency. 
Such  avoidance  also  facilitates  software  modularity,  i.e.,  the  discrete-element  analyzers  for 
the  various  systems  involved  need  not  be  tightly  coupled,  as  the  time  integration  may 
proceed  by  means  of  simple  response-vector  transfers. 

Unfortunately,  these  benefits  exact  a  price:  “conditional  stability"  with  respect  to 
time  increment,  i.e.,  in  order  to  avoid  numerically  unstable  computations,  the  size  of  the 
increment  between  solution  steps  must  be  smaller  than  a  value  determined  by  the  highest- 
frequency  component  of  the  global  equations.  For  “stiff  coupled  systems”,  i.e.,  systems 
with  widely  disparate  frequency  characteristics,  this  drawback  is  truly  burdensome.  In 
fact,  extrapolation  algorithms  can  often  produce  an  integration  method  that  is  uncondi¬ 
tionally  unstable,  i.e.,  one  for  which  no  time  increment  is  small  enough  to  prevent  numerical 
instablity. 

Global  Implicit  Integration:  This  second  method  can  remove  the  burden  of  conditional 
stability  while  preserving  a  relatively  straightforward  approach.  The  main  price  it  exacts 
for  this  is  computational  efficiency,  requiring  the  simultaneous  solution  of  many  algebraic 
equations  poorly  structured  for  rapid  processing.  In  fact,  the  global  equations  are  often 
structured  so  as  to  preclude  practical  computation  by  implicit  algorithms.  In  addition,  the 
requirement  of  global  implicit  integration  for  simultaneous  equation  solution  undermines 
software  modularity. 

Modal  Equation  Integration :  The  method  of  modal  decomposition  overcomes  the  problem 
of  computational  efficiency  by  expressing  the  response  of  the  discrete-element  model  in 
terms  of  “modes"  governed  by  modal  equations  of  motion  whose  number  is  much  smaller 
than  the  number  of  global  equations.  Either  explicit  or  implicit  numerical  time  integra¬ 
tion  may  be  used  to  solve  the  modal  equations,  for  neither  small  time  increments  nor 
simultaneous  equation  solution  is  burdensome  for  small  problems. 

The  price  exacted  by  modal  decomposition,  however,  is  a  second  modelling  stage 
that  is  often  characterized  by  uncertainties  more  troublesome  than  those  associated  with 
the  discrete-element  modelling  stage.  For  example,  while  satisfactory  modes  for  a  linear 
structure  may  often  be  selected  from  a  set  of  eigenvectors  for  an  associated  free-vibration 
problem,  satisfactory  modes  (i.e.,  generalized  functions)  for  a  nonlinear  structure  are  hard 
to  divine. 

Partitioned  Analysis:  This  last  method  achieves  both  computational  efficiency  and  software 
modularity  by  grouping  the  global  equations  into  logical  sets  according  to  computational 
attributes.  A  solution  scheme  is  then  devised  to  optimize  computational  efficiency  subject 
to  stability  and  accuracy  requirements.  The  price  exacted  to  secure  these  benefits  is 
complexity  of  implementation.  The  formulation  and  implementation  of  solution  schemes 
that  are  both  stable  and  accurate  require  techniques  that  are  currently  in  an  embryonic 
stage  of  development,  as  discussed  in  succeeding  sections.  Partitioned  analysis  also  requires 


new  methods  of  software  implementation,  as  discussed  next. 

1.3  SOFTWARE  IMPLEMENTATION. 

Most  codes  used  to  treat  single-system  mechanics  problems  of  engineering  significance 
consist  of  thousands  of  executable  lines  of  Fortran.  Each  is  usually  the  product  of  a 
single  individual  or  of  a  tightly  knit  group  of  developers.  The  treatment  of  coupled-system 
problems  at  the  same  level  of  complexity  may  be  accomplished  in  two  ways: 

1.  Extension  of  an  existing  code  for  one  set  of  field  equations  to  accomodate  discretization 
of  the  additional  field  equations  and  solution  of  these  with  those  of  the  original  code; 
this  is  the  monolithic  software  approach. 

2.  Integration  of  separate  codes  through  flexible  runstream  control  and  data  manage¬ 
ment;  this  'is  the  integrated  software  approach  [8]. 

In  cases  where  one  set  of  field  equations  is  quite  complex,  and  where  the  other  field  equa¬ 
tions  are  quite  simple,  such  as  those  for  viscous  dampers,  the  monolithic  software  approach 
is  to  be  preferred.  However,  when  two  or  more  of  the  field-equation  sets  are  complex,  the 
prospect  of  embedding  the  large  code  for  one  set  into  that  for  another  set  is  awesome. 
Furthermore,  the  effort  will  probably  fail,  for  the  level  of  complexity  escalates  to  a  point 
where  mere  mortals  cam  no  longer  cope  with  it. 

Even  though  the  integrated  software  approach  largely  circumvents  the  horrors  of  monolithic 
code  development,  it  offers  stiff  challenges  of  its  own.  The  integration  of  disparate  soft¬ 
ware  through  the  use  of  runstream-control  and  data-management  utilities  is  in  its  infancy, 
strained  by  the  conflicting  demands  of  software  modularity  and  computational  efficiency. 
Nevertheless,  it  is  apparently  the  only  practical  way. 

x.4  SPSS  PROBLEMS. 

SPSS  is  tasked  to  address  some  of  the  most  difficult  technical  problems  of  modern  tech¬ 
nology.  These  problems  often  involve  extreme  environments  (large  pressures,  high  tem¬ 
peratures),  complex  behavior  (large  motions,  nonlinear  material  response),  and  coupled 
systems  (submerged  or  embedded  structures).  The  people  and  codes  required  to  analyze 
these  problems  are  highly  specialized,  constituting  islands  of  expertise  that  tend  to  resist 
integration.  Partitioned  analysis,  as  a  method  of  computation,  and  integrated  software,  as 
a  method  of  implementation,  provide  perhaps  the  least  painful  means  to  bring  the  people 
and  their  codes  together. 


SECTION  2 

TOOLS  OF  PARTITIONED  ANALYSIS 


2.1  TEMPORAL  DISCRETIZATION. 

In  this  report  we  focus  attention  on  coupled  systems  governed  by  second-order  matrix 
differential  equations,  as  these  are  dominant  in  mechanical  applications.  The  global  form 
of  these  equations  is: 


Mu  +  Cu  +  Ku  =  f  (1) 

where  the  global  solution  vector  u  embodies  all  degrees  of  freedom  of  the  coupled  system, 
M,  C  and  K  are  analogues  of  the  mass,  damping  and  stiffness  matrix,  respectively,  of 
structural  analysis,  and  superposed  dots  denote  temporal  differentiation.  Any  nonlinear 
term  is  assumed  to  be  collected  in  the  forcing  vector  f. 

Transformation  to  First  Order 

A  systematic  treatment  of  the  global  equations  is  considerably  simplified  by  transforming 
the  second-order  system  (l)  into  a  set  of  first-order  equations.  To  this  effect,  and  following 
Refs.  [5,14],  introduce  the  general  auxiliary  vector  v: 

v  =  AMu  +  Bu  (2) 

where  A  and  B  are  square  matrices  (arbitrary  except  for  A  being  nonsingular)  that  may 
be  chosen  so  as  to  simplify  the  implementation  or  to  minimize  the  computational  cost. 
Differentiating  (2)  once  and  substituting  for  u  from  (l)  yields 


v  =  A  (f  -  Ku  -  Cu)  +  Bu 


(3) 


Equations  (2)  and  (3)  can  be  expressed  as  a  set  of  first  order  matrix  equations: 
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where  I  and  0  denote  the  identity  and  null  matrix,  respectively.  To  proceed  with  the 
temporal  discretization,  we  have  to  introduce  a  pair  of  integration  formulas,  one  for  u  and 
one  for  v. 


Time  Integration  Formulas 

We  restrict  our  attention  to  first-order  linear  multistep  formulas  of  the  type 

m  m 

Un  +  ^2  a3Un-3  =  ^  ^  n~ J 

,=1  i=1  (5) 

m  m  v  ' 

vn  +  aiV»-J  =  ^Vn  +  h  03  V»-i 

J=1  J=1 

where  h  is  the  integration  stepsize,  6  =  (3q}i  is  a  scaled  stepsize,  ay  and  /?y  are  integrator 
coefficients,  n  is  the  current  time  step  index,  and  m  is  the  number  of  steps  used  in  the 
construction  of  the  formulas.  Note  that  if  <5  =  0,  then  u„  and  vn  can  be  directly  computed 
in  terms  of  known  past  quantities.  Thus,  formulas  in  which  thecoefficient  /?0  vanishes  are 
called  explicit,  otherwise  they  are  called  implicit.  Various  generalizations  of  (5)  are  possible 
but  not  treated  here.  * 

Difference  Equations 

Introduction  of  (5)  into  (4)  yields  the  difference  equations 
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in  which  the  historical  vectors  h„  are  given  by 
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K  =  £(-ayU„_y  +  h(3jiin-j)  (7) 

3=1 

m 

K  =  ^2{~a3^n-3  +  ^Vn_y)  (8) 

3=1 

If  6  ^  0,  the  auxiliary  vector  v  can  be  eliminated  from  (7)  to  produce  the  implicit  system 

Eun  =  g„  (9) 


*  Two  immediate  extensions  are  possible,  (l)  Different  formulas  may  be  used  for  u  and  v  as 
proposed  in  [5];  this  extension  has  not  proven  useful,  however.  (  Different  formulas  may  be 
used  for  different  components  of  the  the  solution  vector  so  that  the  coefficients  a,  and 
vary  across  the  system;  this  extension  has  been  advantageously  used  in  fluid-structure  inter¬ 
action  analysis.  Introduction  of  these  more  general  forms  is  straightforward,  but  complicates 
subsequent  expressions. 
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where 

E  =  M  +  6C  +  S2K  (10) 

gn  =  [M  +  (C-A-1B)]h“  +  5Ah:  +  ^fn  •  (11) 

Equation  (9)  is  a  linear  algebraic  system  that  contains  the  same  number  of  equations  as 
the  global  coupled  system  (1).  This  system  must  be  solved  at  each  time  step.  To  complete 
the  formulation,  however,  the  two  arbitrary  matrices  A  and  B  must  be  selected. 

Selection  of  Auxiliary  Matrices 

Two  widely  used  choices  for  A  and  B  in  (2)  are 

v=u  (A  =  M~ l,  B  =  0)  (12) 

v  =  Mu  +  Cu  (A  =  I,  B  =  C)  (13) 

In  [5-7]  these  choices  are  designated  as  the  conventional  and  Jensen’s  [14]  choices,  respec¬ 
tively.  They  result  in  different  appearances  of  the  difference  equations.  For  example,  the 
choice  (13)  reduces  E  and  gn  in  (10-11)  to 


! 


E  =  M  +  6C  +  £2K 
gn  =  (M+*c)h:+«i:+*2fn 


2.2  PARTITIONING. 


Motivations  for  undertaking  partitioned  analysis  methods  are  discussed  in  Section  1.2. 
This  section  reviews  the  mechanics  of  the  process.  Partitioning  involves  two  operations: 
splitting  and  blocking. 

Splitting 

Partitioned  analysis  of  second  order  systems  relies  on  the  following  additive  decomposition 
of  the  damping  and  stiffness  matrices:  * 

D  =  D7  +  De 

/  B  (16) 

k  =  k7  +  k£: 

where  superscripts  /  and  E  stand  for  “implicit  part”  and  “explicit  part”,  respectively,  on 
account  of  the  interpretation  given  later.  The  induced  decomposition  of  the  coefficient 
matrix  E  of  the  implicit  difference  equation  (10)  is  ** 

*  The  limiting  cases  of  globally  explicit  and  globally  implicit  treatment  correspond  to  KJ  = 
Dr  =  0  and  KE  =  Ds  =  0,  respectively. 

**  The  procedure  described  herein  is  known  as  algebraic  partitioning  and  is  described  in  further 
detail  in  Refs.  [7,19,20],  There  is  an  alternative  procedure,  called  differential  partitioning, 
in  which  partitioning  is  carried  out  at  the  differential  equation  level  with  interaction  terms 
viewed  as  predicted  pseudo-forces.  This  was  historically  the  first  procedure  to  be  used  for 
fluid-structure  interaction  [17],  but  is  presently  less  important  than  algebraic  partitioning. 


raw 


-  •  .  'j  V  V  V  V  V//.  ✓  .  (  •r,  »\  r.. 


E  =  E  +  E 


(17) 


in  which 


EJ  =  M  +  6DJ  +  SK1 
Ee  =  6T>e  +  6Ke 


(Note  that  the  mass  matrix  M  is  not  decomposed.)  The  term  involving  Ee  is  transferred 
to  the  right-hand  side,  and  a  predictor  is  applied: 

EJun  =  g„  —  EM  (19) 

where  gn  is  given  by  (11),  and  u£  is  a  value  predicted  from  previous  solutions.  A  general 
form  of  this  predictor  is 

m 

=  ^(dtn-jUn-y  +  6$n-jXln-j  +  62*)n-j-  Un_y)  (20) 

y=i 

where  the  derivative  terms,  if  present,  are  those  calculated  in  the  solution  process. 

To  complete  the  partitioning  scheme,  the  block  structure  of  matrices  and 

Ke  has  to  be  specified  as  suggested  by  the  physics  of  the  problem  as  well  as  computer 
implementation  considerations.  This  topic  is  treated  below  in  an  example-oriented  form, 
starting  from  the  simplest  case  of  two  coupled  fields. 

Two-field  Problems 

The  simplest  coupled  system  involves  two  interacting  fields  (AT,  Y).  Partition  u  as 

“={«;}  <2i> 

where  vector  ux  contains  the  degrees  of  freedom  of  field  X  and  vector  uy  the  degrees  of 
freedom  of  field  Y .  The  corresponding  form  of  (1)  is  * 

Mx  0  f  ux  1  I"  T)xx  Djy  j  ux  1  r  KZI  Kjy  1  (  ux  \  _  f  fx  1  (Or)\ 

0  MyJ\uJ  +  [Dyx  DyyJ  \uy/+  [Kyx  KyyJ\uy/-\fy/  ^ 

This  type  of  coupled  system  arises  naturally  in  contexts  in  which  the  state  variables  ux  of 
field  X  are  conjugate  (in  the  virtual  work  sense)  to  those  of  field  Y.  ** 


*  Note  the  absence  of  acceleration  coupling  terms;  this  is  essential  for  the  success  of  these 
methods.  It  is  shown  in  [7]  that  such  terms  can  always  be  eliminated  should  they  be  present 
in  the  original  problem. 

**  This  happy  circumstance  occurs  naturally  in  the  case  of  a  structure  interacting  with  an 
external  medium.  For  example,  if  the  medium  is  an  acoustic  fluid,  u*  and  u,  become 
structural  nodal  velocities  and  fluid  surface  pressures,  respectively. 
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Even  for  only  two  fields,  there  is  a  surprisingly  large  number  of  different  block  parti¬ 
tions:  ten  to  be  exact.  The  possible  configurations  of  KE  are: 


0  0 
0  0 


o  Kxy 

Kyx  0 


0 

0 


0 

Kyx 


0 

0  ol 

r  o 

1  0 

0 

Kxy' 

Kyy. 

O 

H 

a 

w 

[KyX  Kyy\ 

0 

Kyy 

Kxy' 

\KXX  1 

»  1 

[Kxx  1 

0  1 

K« 

Kxy 

Kyy. 

L  0  K-yy  J 

[Kyx  KyyJ 

Kyx 

KyyJ 

(23) 


with  a  similar  partition  assumed  for  DB.  (There  are  actually  4  2  =  16  block  partitions, 
but  six  correspond  to  trivial  switchings  of  X  and  Y.) 


Terminology 

Once  a  two-field  partition  is  selected,  field  X  is  said  to  be  treated  explicitly  if  the  diagonal 
block  Kxx  ends  up  in  KB,  and  implicitly  otherwise.  Similarly,  field  Y  is  said  to  be  treated 
explicitly  if  the  diagonal  block  Kyy  ends  up  in  KE ,  and  implicitly  otherwise. 

A  partition  in  which  at  least  one  field  is  treated  explicitly  is  called  implicit-explicit 
(I-E).  If  all  fields  are  treated  implicitly,  the  partition  is  called  implicit-implicit  (I-I).  In  an 
I-E  partition,  the  explicit  part  is  determined  simply  by  extrapolation  whereas  in  an  I-I 
partition  a  final  implicit  pass  is  required.  An  I-I  partition  is  easily  recognized  by  the  lack 
of  diagonal  blocks  in  KE  and  DE.  (This  terminology  holds  for  any  number  of  fields.) 

In  the  case  of  the  two-field  partitions  (23),  the  first  and  last  correspond  to  the  limit 
cases  of  globally  implicit  and  globally  explicit  treatment  of  the  problem,  respectively.  Of 
the  other  eight,  two  (the  third  and  sixth  one)  are  implicit-implicit,  one  (the  eighth  one) 
is  explicit-explicit  and  thus  trivially  equivalent  to  the  last  one;  and  the  other  five  are 
implicit-explicit. 

The  most  practically  important  is  the  third  one,  the  staggered  I-I  partition,  which 
has  been  extensively  used  in  the  structure-fluid  interaction  problems  (Section  4.1).  The 
most  interesting  I-E  partition  is  the  second  one,  which  follows  under  the  purview  of  the 
“DOF-by-DOF”  partition  examined  later. 


Simplified  Three-Field  Problem 

When  more  them  two  fields  are  involved,  the  number  of  possible  partitions  increases  very 
rapidly.  To  give  the  flavor  of  this  “combinatorial  explosion”  consider  a  problem  in  which  we 
have  two  physical  fields:  X  and  Y,  which  interact  through  a  boundary  B.  The  boundary 
is  viewed  as  a  separate  field,  so  the  partition  of  the  solution  vector  is 


u  = 


Ul 
{  u«> 


(24) 


Physically  this  case  corresponds  to  the  interaction  of  two  domain-type  discretizations,  e.g. 
two  finite  element  meshes.  For  such  problems  the  boundary  unknowns  warrant  special 
treatment. 
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The  K  matrix  for  the  coupled  system  has  the  block  structure 


Kzz  K  xb  0 

Kb*  K^  +  K*,  Kby 

0  Kb*  Kyy 


(25) 


with  a  similar  structure  assumed  for  matrix  D. 

There  are  102  distinct  partitions  of  K,  most  of  which  are  of  little  interest.  In  what 
follows  we  examine  a  few  that  have  proven  practically  useful. 


Node-by-Node  I-E  Partition 
Consider  the  partition 


K** 

K*b 

0 

0 

0 

0 

K  = 

Kb* 

K£,  +  K&, 

Kby 

+ 

0 

0 

0 

0 

0 

0 

0 

Kyb 

Kyy. 

(26) 


The  state  vector  for  Y  is  computed  explicitly;  then  the  boundary  and  X-field  unknowns 
are  computed  by  a  fully-implicit  scheme. 

This  I-E  partition  was  proposed  by  Belytschko  and  Mullen  [1-3],  who  described  it  in 
a  more  physical  context.  It  is  a  natural  partition  for  two  interacting  finite  element  meshes 
if  boundary  nodes  are  identified  as  part  of  the  implicitly-treated  mesh.  This  interpretation 
motivates  its  name. 


Element-by-Element  I-E  Partition 


The  node-by-node  partition  views  boundary  unknowns  as  part  of  the  implicit  X  field;  the 
resulting  "partition  anisotropy"  is  reflected  in  an  unsymmetric  KE .  Now  consider  the 
symmetric  partition 


K  = 

Kx* 

Kb* 

K*b 

K£> 

0 

Kby 

+ 

1 

O  O 

0 

Kbb 

0 

0 

0 

0 

0 

L° 

Kyb 

Kyy  _ 

(27) 


Field  X  and  its  connection  to  the  boundary  is  treated  implicitly  whereas  field  Y  and  its 
connection  to  the  boundary  are  treated  explicitly.  Thus,  part  of  the  boundary  values  are 
obtained  explicitly,  and  partly  implicitly. 

This  partition  is  natural  for  finite  element  codes  in  which  it  is  desired  to  label  elements 
as  implicit  or  explicit  as  proposed  by  Hughes  and  Liu  [12-13].  The  fact  that  boundary  nodes 
need  not  be  explicitly  tagged  as  such  simplifies  programming.  On  the  other  hand,  if  the 
interface  between  the  X  and  Y  fields  is  extensive,  this  partition  entails  more  computational 
work  than  a  node-by-node  partition. 
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DOF-by-DOF  I-E  Partition 

As  final  I-E  example,  consider  the  symmetric  partition 


K** 

Kxb 

0  1 

0 

0 

Kfcy  + 

0 

0 

0 

Kyb 

0 

J 

0 

0 

Here  only  the  V-field  vector  is  treated  in  explicit  fashion;  note  the  difference  with  regard 
to  the  node-by-node  partition. 

This  partition,  proposed  in  [19],  is  useful  when  different  degrees  of  freedom  are  to  be 
treated  by  different  algorithms.  For  example,  rotational  degrees  of  freedom  associated  with 
beam  or  shell  structures  display  higher  frequencies  than  translational  freedoms.  It  would 
be  advantageous  to  treat  the  former  implicitly  whereas  the  latter  are  treated  explicitly. 

Staggered  /-/  Partition 

Although  the  staggered  partition  has  been  primarily  used  in  two-field  coupled  problems, 
it  can  also  be  formulated  for  the  two-field-plus-boundary  problem  as 

K «  Kxb  0  ]  [0  0  O' 

K=  Kj*  K£,  K*  +  0  0  0  (29) 

0  0  Kyy  ]  [o  Kyi>  0 

Since  all  diagonal  blocks  remain  on  the  left  side,  the  partition  is  implicit-implicit.  Boundary 
values  have  to  be  predicted  for  field  Y ,  but  the  determination  of  uy  relies  on  a  implicit 
(corrector)  pass. 

2.3  CHARACTERISTIC  EQUATION. 

The  characteristic  equation  determines  stability  and  accuracy  characteristics  of  a  linearized 
partitioned  analysis  procedure.  It  is  presented  here  (without  proof)  as  a  backup  to  some 
statements  made  in  Section  2.4.  The  derivation,  which  is  covered  in  further  detail  in 
[7,19,21],  begins  with  the  standard  assumption 

un  =  Au„_x  (30) 

where  A  is  a  complex  number  called  the  amplification  factor.  Denote  the  characteristic 
polynomials  of  the  integration  formulas  (5)  and  predictor  (22)  by 

m 

P(A)  =  A™  +  £«/Am'y 
/=! 
m 

a(A)  =  Am  +  £/?yA— i  (3i) 

y=i 

171  m  771 

«W  = 

y=i  y=  i  y=i 
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With  these  definitions,  the  characteristic  equation  may  be  expressed  as 


C(A)  =  Cu{\)  +  CP{\)  (32) 

where 

Cu{  A)  =  p2  M  +  SpcT>  +  62o2K 

Cp( A)  =  (e  -  \m)(6cDDE  +  62ck Ke)  ^ 

in  which  cD  and  cK  are  functions  of  p  and  a  that  depend  on  the  computational  path 
(Section  2.4).  Cu { A)  is  the  characteristic  equation  for  a  global  (unpartitioned)  treatment  of 
the  global  system,  while  Cp( A)  is  a  correction  term  that  brings  the  effect  of  the  partitioning 
process. 

2.4  IMPLEMENTATION. 

Once  the  partition  is  selected,  three  more  ingredients  must  be  selected  to  complete  the 
algorithm  and  its  computer  implementation: 

•  An  specific  integration  formula 

•  A  predictor  formula 

•  The  computational  path 

These  three  aspects  are  not  independent.  Much  of  the  progress  made  to  date  in  partitioned 
analysis  theory  consists  of  a  the  study  of  these  interrelations  for  certain  classes  of  problems 
of  interest  in  computational  mechanics.  Any  progress  in  this  direction  helps  to  alleviate 
the  forbidding  combinatorial  nature  of  the  design  process.  The  following  material  reflects 
our  present  understanding  and  may  change  or  expand  in  the  future  as  the  theory  evolves. 

Selecting  the  Integrator 

If  the  global  system  is  linear  and  undamped  or  lightly  Rayleigh-damped  the  trapezoidal 
rule  (TR) 

un  -  u„_ !  =  {h/ 2) (un  +  un-t) 

Vn-yn-i  =  (h/2)(vn  +  vn-t) 

is  strongly  recommended.  The  TR  has  the  smallest  truncation  error  among  all  A-stable* 
formulas  of  the  linear  multistep  type,  introduces  no  numerical  damping,  and  being  a  one- 
step  method  it  is  self-starting.  For  explicitly  treated  systems  in  I-E  partitions,  (34)  becomes 
the  central  difference  method,  which  is  another  old  favorite  for  linear  undamped  systems. 
The  main  proviso  is  to  avoid  certain  computational  paths  that  cause  the  TR  to  generate 
spurious  “beating”  oscillations  [18]. 

As  we  depart  from  the  “undamped  oscillator”  model,  the  situation  becomes  progres¬ 
sively  less  clear.  If  damping  mechanisms  of  radiation  type  enter  the  picture  the  coupled 


*  An  integration  formula  is  called  A-stable  if  it  is  stable  for  any  stepsize  when  applied  to  the 
test  equation  u  =  Au  in  which  R(A)  >  0. 
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system  exhibits  more  parabolic  characteristics,  and  the  TR  may  not  be  necessarily  the 
best  scheme,  or  even  a  stable  one. 

If  nonlinearities  are  present,  it  is  well  known  that  the  TR  may  fell  prey  to  spurious 
oscillations.  If  the  nonlinearities  are  mild  in  nature,  changing  the  TR  to  a  closely  related 
formula,  the  midpoint  rule,  may  be  sufficient  to  eliminate  the  oscillations.  For  strong 
nonlinearities,  however,  it  is  generally  needed  to  introduce  artificial  damping  to  stabilize 
the  integration.  This  damping  may  be  introduced  in  the  integrator  formula**,  or  added 
to  the  governing  differential  equations  as  a  stepsize-dependent  dissipative  term*** 

In  summary:  for  generally  damped  and/or  strongly  nonlinear  systems  a  comprehen¬ 
sive  theory  is  lacking,  and  we  must  resort  to  computer-aided  experimentation  to  design 
partitioned  analysis  procedures. 

Selecting  the  Predictor 

The  next  task  pertains  to  to  the  selection  of  u£  in  (19).  A  body  of  theory  as  regards 
predictors  of  the  linear  multistep  type  (20)  presently  exists. 

In  the  initial  development  of  partitioned  analysis,  predictors  were  restricted  to  com¬ 
binations  of  past  solutions  (terms  with  the  a}  coefficients).  Computed  derivatives  were 
added  later  to  compensate  for  the  effect  of  the  computational  path,  and  to  attain  optimal 
accuracy  in  representation  of  certain  motions  [20] . 

Computational  Path 

Four  computational  paths,  labeled  (O’),  (0),  (l)  and  (,  may  be  followed  in  advancing 
the  solution  over  a  typical  time  step.  The  general  expression  of  such  paths  is  flow-charted, 
e.g.  in  [5-7].  The  essential  difference  between  the  paths  is  the  manner  in  which  the  auxiliary 
vector  v„  and  its  temporal  derivative  vn  are  calculated  at  each  step. 

The  path  identification  index  (0,  1  or  2)  gives  the  number  of  backward-difference 
operations  performed  in  the  determination  of  u„,  v„  and  v„  in  each  time  step.  In  Ref.  [19] 
it  is  shown  that  this  index  plays  an  important  role  in  the  computational  error  propagation 
behavior. 

If  the  global  system  is  not  treated  by  a  partitioned  analysis  method,  the  choice  of 
computational  path  affects  only  computational  aspects:  implementation  efficiency  and 
error  propagation.  Stability  and  algorithmic  accuracy  are  not  affected. 

But  if  a  partitioned  analysis  procedure  is  adopted,  the  characteristic  equation  is  found 
to  depend  on  the  computational  path  as  the  expression  (33)  makes  evident.  This  was  a 
key  research  result.  The  practical  consequence  is  that  a  carefully  stabilized  implemen¬ 
tation  may  suddenly  become  unstable  if  apparently  innocuous  changes  are  made  to  the 
computational  sequences.  This  is  a  grave  concern  since  application  software  is  frequently 
updated  to  meet  new  problems. 

**  For  example,  the  3-step  Park  method  [16],  or  the  2-step  Gear  method  [14].  This  route  is 
commonly  followed  in  implicitly  treated  fields. 

***  This  is  normally  done  for  explicitly  treated  fields,  for  example  the  fluid  volume  system  in  the 
cavitation  analysis  discussed  in  Section  4.3. 


Path  Compensation 

The  dependence  of  the  characteristic  equation  (32)  on  the  computational  path  is  through 
terms  cD  and  cK.  The  dependence  on  the  predictor  is  through  the  term  e.  The  particular 
form  of  this  dependence  made  two  key  results  possible:  (1)  it  is  possible  to  adjust  predictors 
so  that  the  same  characteristic  equation  results  for  all  computational  paths,  and  (  these 
adjusted  predictors  are  implementable ,  that  is,  of  the  form  (20),  if  some  mild  constraints 
are  enforced.  This  key  result  is  further  discussed  in  [7]. 

These  results  have  important  practical  consequences.  If  a  partitioned  analysis  pro¬ 
cedure  is  stabilized  for  one  computational  path,  it  can  be  stabilized  for  others  by  simply 
changing  the  predictor  formula  as  given  by  theory. 
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SECTION  3 

DESIGN  TECHNIQUES 


3.1  DESIGN  CRITERIA. 

The  design  of  a  partitioned  analysis  procedure  aims  to  find  the  optimal  combination  of  gov¬ 
erning  equations,  partitioning  scheme,  time  integration  algorithm,  extrapolation  algorithm 
and  computational  path  that  conforms  to  given  constraints. 

The  adjective  “optimal”  can  be  characterized  as  follows.  The  optimal  procedure 

enjoys  unconditional  stability; 
is  highly  accurate; 

is  easy  to  implement  on  the  computer; 

is  economical  to  run;  and 

satisfies  software  modularity  constraints. 

Unfortunately,  the  stated  goals  are  mutually  contradictory.  Instead  of  searching  for  the 
Holy  Grail  we  have  to  settle  on  a  compromise.  The  nature  of  the  compromise  is  dictated 
by  the  relative  weight  attached  to  the  procedure  rating  attributes  (stability,  accuracy, 
simplicity,  efficiency)  as  well  as  the  modularity  constraints. 

Relative  attribute  weighting  is  seldom  stated  a  priori,  but  comes  implicitly  from 
answering  the  question:  which  use  will  the  analysis  procedure  be  put  to?  Some  examples:  if 
tracing  a  late-time  or  steady-steate  response  is  deemed  important,  attaining  unconditional 
stability  should  be  given  high  priority.  If  the  resulting  computer  program  is  to  be  heavily 
used  in  a  production  environment,  or  complex  nonlinear  problems  attacked,  computational 
efficiency  might  be  stressed  over  ease  of  implementation.  High  accuracy  is  rarely  crucial 
when  solving  problems  in  structural  mechanics,  but  may  be  required  in  exotic  applications 
such  as  orbit  calculations  or  aiming  of  electromagnetic  beams. 

The  foregoing  design  criteria  are  next  examined  in  more  detail. 

3.2  MODULARITY. 

A  common  design  constraint  is  related  to  software  modularity:  the  search  for  an  optimal 
partitioned  analysis  procedure  has  to  account  for  the  prior  existence  of  one  or  more  com¬ 
puter  programs  that  model  and  solve  part  of  the  complete  problem.  The  idea  is  to  make  use 
of  such  programs  as  components  of  the  software  system  that  solves  the  coupled  problem. 
Why?  As  discussed  in  Section  1,  existing  codes  may  represent  substantial  investments  in 
technology  development  and  user  training.  One  of  the  strength  of  partitioned  analysis  is 
the  relative  ease  with  which  solution  procedures  can  be  adapted  to  meet  these  modularity 
requirements. 

The  modularity  constraint  impacts  primarily  the  selection  of  the  governing  coupled 
equations  and  the  freedom  with  which  these  equations  can  be  manipulated  to  improve 
numerical  stability  characteristics.  On  a  secondary  level,  it  may  curtail  the  selection  of 
integration  algorithm,  extrapolation  and  computational  path. 


3.3  STABILITY. 


Numerical  stability  is  the  toughest  algorithmic  property  to  control  in  partitioned  analysis, 
and  success  or  failure  hinges  heavily  on  attaining  satisfactory  stability  characteristics. 

But  what’s  “satisfactory”?  We  would  obviously  like  to  achieve  unconditional  stability , 
i.e.  for  any  time  stepsize.  But  this  statement  is  equivocal,  as  shown  by  the  following 
counterexamples. 

Conditionally  Stable  Subsystem.  Suppose  it  is  known  a  priori  that  a  system  component 
has  to  be  treated  by  an  explicit  integration  algorithm  that  limits  the  maximum  stepsize 
to  hmax.  This  situation  may  arise  because  of  software  modularity  requirements,  faster 
implementation  under  schedule  constraints,  or  simply  the  computational  impracticality 
of  using  an  implicit  integrator.  Now  the  stable  stepsize  for  the  coupled  system  cannot 
exceed  hmax,  for  the  stability  of  the  whole  cannot  exceed  that  of  the  (isolated)  parts.  It 
is  therefore  pointless  to  ask  for  unconditional  stability.  The  correct  objective  is  to  avoid  a 
reduction  of  the  stability  limit  when  the  coupled  system  is  treated. 

Physical  Instability.  If  the  application  intends  to  uncover  physical  instability,  for  example 
flutter  in  aeroelastic  analysis,  numerical  “overstability”  is  undesirable.  More  specifically, 
many  highly  stable  time  integration  algorithms  display  heavy  numeric  damping.  This 
damping  may  mask  or  even  suppress  the  onset  of  physical  stability;  so  it’s  possible  to  have 
too  much  of  a  good  thing.  If  the  detection  of  physical  instability  is  deemed  important,  an 
additional  objective  is  to  have  the  boundaries  of  numerical  and  physical  instability  coincide 
independently  of  stepsize.  If  this  objective  cannot  be  met,  there  should  be  provisions  for 
iterating  on  the  coupled  system  while  keeping  the  time  frozen;  physical  instability  is  then 
revealed  by  lack  of  convergence. 

In  light  of  the  foregoing  points,  the  design-for-stability  should  try  to  (l)  minimize  the 
degradation  of  numerical  stability  with  respect  to  the  uncoupled-system  stability  limit, 
and  (2)  minimize  masking  of  physical  dynamical  instability  if  the  latter  is  of  concern. 

3.4  ACCURACY. 

Accuracy  characterizes  the  fidelity  with  which  the  computed  solution  approximates  the 
exact  solution  of  the  semi-discrete  governing  equations.  We  are  therefore  talking  about 
accuracy  in  the  time  domain,  and  not  about  errors  introduced  by  the  discretization  in 
space  that  produces  a  system  of  finite  number  of  degrees  of  freedom. 

The  distinction  helps  in  explaining  why  very  high  accuracy  requirements  are  uncom¬ 
mon  in  computational  mechanics.  For  problems  of  engineering  significance,  a  global  ac¬ 
curacy  of  1%  to  10%  in  peak  displacements  and  velocities  is  generally  satisfactory,  since 
uncertainties  of  this  order  axe  inherent  in  the  spatial  discretization,  material  modelling 
(especially  constitutively  nonlinear  behavior)  and  representation  of  external  forces. 

Accuracy  under  1%  is  only  required  in  simple  test  problems  used  to  verify  the  method 
and  its  implementation,  and  in  applications  where  certain  “system  integrated”  character¬ 
istics  (motion  of  center  of  mass,  acoustic  cross  section,  etc.)  are  of  primary  interest. 

It  follows  that  accuracy  considerations  are  not  a  driving  life-and-death  proposition 
in  method  design,  the  way  stability  is.  The  proper  time  to  consider  it  is  when  several 
competing  stable  formulations  remain  after  modularity  constraints,  if  any,  are  met. 


The  final  product  of  a  thorough  accuracy  analysis  should  be  a  comprehensive  set 
of  rules  for  controlling  the  time  stepsize,  perhaps  automatically.  This  is  seldom  done  in 
practice,  for  heuristic  stepsize  selection  rules  based  on  the  physics  of  the  problem  work 
surprisingly  well.  It  has  been  our  experience  that  for  shock-excited  structural  systems 
of  engineering  interest  the  natural  stepsize  required  to  capture  the  characteristics  of  the 
excitation  yield  early-time  responses  of  comfortable  accuracy.  This  is  of  course  typical  of 
integration  processes  in  spatially  elliptic  systems,  which  smooth  the  input  data.  Intuition 
may  fail,  however,  if  the  system  under  consideration  is  not  of  elliptic  type. 

3.5  EASE  OF  IMPLEMENTATION. 

This  attribute  relates  to  the  degree  of  difficulty  with  which  a  partitioned  analysis  procedure 
can  be  translated  into  a  working  computer  program.  This  is  not  an  easily  quantifiable 
attribute,  but  rather  a  set  of  guidelines  that  should  be  kept  in  mind  throughout  the  design 
process. 

If  modularity  constraints  apply,  implementation  difficulty  raises  according  to  the 
amount  and  nature  of  modifications  that  must  be  made  to  the  existing  program(s),  as 
itemized  below. 

Bare-Minimum  Modifications.  Make  provisions  to  receive  input  interaction  data  and  return 
output  interaction  data  at  each  time  step.  This  is  straightforward  if  all  interaction  data  can 
be  organized  and  transmitted  in  vector  or  diagonal-matrix  form;  more  difficult  if  general 
matrix  structures  are  involved. 

Moderate  Modifications.  This  case  is  typified  if  the  procedure  implementation  calls  for 
innocuous  interacting-field  terms,  for  example  a  diagonal  damping  matrix,  to  enter  the 
left-hand-side  of  the  programmed  equations.  Or  if  additional  routines  must  be  to  written 
to  calculate  derived  quantities,  such  as  a  residual  dynamic  force  vector,  that  must  be 
passed  to  other  codes.  The  key  point  is  that  the  basic  solution-advancing  algorithm  is  not 
affected. 

Substantial  Modifications.  This  is  typified  by  an  implementation  that  requires  the  intro¬ 
duction  of  a  consistent  mass,  non-diagonal  damping  or  unsymmetric  stiffness  matrix  in  a 
structural  analysis  code  not  equipped  to  handle  such  terms.  It  thus  entails  substantial  and 
expensive  modifications  to  innermost  code  levels.  Modifications  of  this  magnitude  should 
be  avoided  if  at  all  possible. 

The  following  general  guidelines  apply  to  new  computer  implementations  that  are  not 
affected  by  modularity  constraints. 

Explicit  vs.  Implicit.  Explicit  time  integration  schemes  are  far  easier  to  implement  than 
implicit  ones.  Furthermore,  automatic  stepsize  control  schemes  are  more  highly  developed 
for  the  explicit  case. 

Treatment  of  Nonlinearities.  If  a  nonlinear  system  is  treated  by  implicit  integration,  im¬ 
plementation  of  a  pseudo-force  scheme  is  easier  than  an  incremental  (sometimes  called 
“tangent  stiffness”)  sc1  me. 

Operator  Symmetry.  Implicit  integration  schemes  that  involve  large  and  sparse  symmetric 
coefficient  matrices  are  easier  to  implement  than  those  that  involve  sparse  unsymmetric 
matrices.  If  the  matrices  are  small  and  dense,  however,  the  difference  is  minor. 


3.6  COMPUTATIONAL  EFFICIENCY. 


The  issue  of  computational  efficiency  can  hardly  be  quantified  with  general  rules.  There  are 
imponderables  such  as  the  host  computer  hardware  and  the  programmer’s  ability,  which 
are  beyond  the  scope  of  this  report. 

When  tackling  a  brand  new  application  of  partitioned  analysis,  it  is  probably  best  not 
to  worry  excessively  about  efficiency.  Chances  are  that  an  initial  implementation  will  be 
primarily  used  by  a  few  people  as  a  research  tool  or  for  analyzing  some  exotic  problems, 
since  it  takes  time  to  build  and  train  a  user  base. 

As  a  program  enters  more  extensive  production  use,  its  computational  efficiency  may 
be  subject  to  closer  scrutiny.  It  is  not  uncommon  then  to  find  that  the  processing  of  a 
single  component  dominates  the  total  computational  cost.  For  extemal-fluid/structure 
interaction  in  which  the  fluid  is  treated  by  a  boundary-element  technique,  the  structural 
analyzer  dominates  the  cost  (especially  if  constitutive  nonlinearities  are  present).  If  the 
structure  is  surrounded  by  a  three-dimensional  fluid-volume  mesh,  the  latter  dominates  the 
overall  cost.  These  observations  may  help  in  the  redesign  of  partitioned  analysis  procedures 
for  subsequent  code  releases. 

3.7  DESIGN  STEPS. 

This  section  describes  stages  of  the  design  of  a  partitioned  analysis  procedure  in  a  step- 
by-step  fashion. 

STEP  o.  Formulation  of  Original  Equations. 

The  selection  of  the  semi-discrete  governing  equations  for  the  coupled  system  is  a  very 
critical  step.  Very  often  these  equations  appear  in  a  mind-boggling  variety  of  forms  based 
on  different  primary  variables,  asymptotic  approximations,  etc.  Choosing  the  right  form 
from  the  ouset  requires  considerable  insight. 

As  a  general  principle,  try  to  select  conjugate  formulations  for  interacting  field  pairs. 
Conjugate  formulations  are  based  on  primary  variables  that  are  conjugate  in  an  energy 
sense,  such  as  velocities  and  pressures. 

If  the  interacting  fields  are  not  based  on  conjugate  variables,  it  may  become  necessary 
to  treat  the  boundary  unknowns  as  a  third  field.  This  can  increase  the  complexity  of 
the  implementation  as  well  as  the  time  spent  in  finding  a  suitable  partitioned  analysis 
procedure. 

If  abundant  software  for  treating  a  particular  subsystem  exists,  the  software  modu¬ 
larity  constraints  discussed  in  Section  3.1  may  help  in  narrowing  the  search  for  governing 
field  equations. 

If  the  problem  contains  nonlinearities,  transfer  them  to  the  right-hand  side  as  a  pseudo- 
force  vector.  The  design  process  cam  only  handle  linear  systems. 

STEP  l.  Temporal  Discretization 

Apply  a  linear  multistep  time  integration  formula  to  the  equivalent  first-order  system  as 
described  in  Section  2.1.  The  end  result  should  be  a  difference  set  of  equations  such  as 
(2.10).  An  experienced  designer  tries  to  keep  the  time  integration  formula  in  symbolic 
form  as  long  as  possible. 


STEP  2.  Partitioning  Scheme 


Partition  the  coupled  difference  system  by  decomposing  K  and  D  as  indicated  in  Sec¬ 
tion  2.2,  and  then  choosing  a  block  structure  dictated  by  the  nature  of  the  problem  or 
modularity  requirements. 

STEP  3.  Select  Integrator,  Extrapolator  and  Computational  Path 

The  theory  of  partitioned  analysis  shows  that  these  three  choices  are  interconnected,  so 
the  freedom  is  not  so  great  as  it  seems.  In  practice  it  often  works  best  to  select  a  simple 
integration  formula,  say  one-step,  which  is  combined  with  the  recommended  extrapolator 
for  a  certain  computational  path,  say  (l).  Results  obtained  by  such  a  choice  can  be  easily 
translated  to  other  combinations  without  having  to  redo  the  stability  analysis. 

STEP  4.  Stability  analysis 

Two  preparatory  tasks  are  involved. 

Reduction  to  Model  Problem.  It  is  possible  to  derive  stability  results  directly  results  directly 
at  the  matrix  level  for  a  limited  class  of  problems,  namely  symmetric,  undamped,  diagonal- 
mass  interacting  systems.  But  generally  an  analysis  at  this  level  is  unwieldy.  A  generalized- 
coordinate  technique  is  applied  to  reduce  the  full  difference  system  reduced  to  a  model 
system  with  as  many  equations  as  fields.  The  stage  at  which  this  reduction  is  performed 
is  not  too  relevant.  Try  to  produce  dimensionless  equations  with  a  minimal  number  of 
independent  parameters. 

Drop  ignorable  terms.  Drop  all  inhomogeneous  terms  from  the  model  system;  this  process 
gets  also  rid  of  nonlinearities.  An  experienced  designer  also  drops  at  this  point  terms 
that  are  known  to  have  no  effect  on  stability;  for  example  structural  damping.  Other 
parameters  may  be  dropped  in  sucessive  passes  of  the  design  process  after  an  initial  pass 
pinpoints  the  critical  parameters. 

Stability  Analysis.  This  ultimately  reduces  to  the  testing  whether  stability  polynomials 
obtained  from  the  characteristic  equations  are  stable.  The  presence  of  problem  parameters 
makes  such  tests  nontrivial.  For  polynomials  of  order  2  or  3  the  tests  may  be  performed  by 
hand  carrying  the  parameters  along.  For  higher  order  the  aid  of  the  computer  is  necessary. 

If  the  analysis  reveals  satisfactory  stability  characteristics  (in  the  sense  of  Section  3.1)  the 
design  process  continues  with  STEP  6. 

If  not,  an  assessment  of  degree  of  stability  is  useful.  This  measures  the  sensitivity  of 
the  stability  boundaries  to  physical  amd  algorithmic  parameters  and  can  be  obtained  by 
perturbing  the  stability  radius  about  unity. 

If  the  boundary  is  sensitive  to  algorithmic  parameters,  for  example  a  predictor  coef¬ 
ficient,  return  to  STEP 3.  If  iterations  produce  no  improvement,  and  there  is  a  choice  of 
partitions,  return  to  STEP  2.  Otherwise  continue  to  the  following  step. 

STEP  5.  Formulation  of  Modified  Field  Equations 

The  governing  field  equations  of  the  same  coupled  problem  can  be  generally  presented  in 
an  infinite  number  of  mathematically  equivalent  formulations.  These  new  formulations  are 


obtained  by  solving  from  terms  from  one  system  and  substituting  into  another  system. 
The  process  is  known  as  augmentation. 

For  problems  such  as  structure/medium  interaction,  the  augmentation  process  can 
radically  change  the  stability  characteristics  of  the  time-integration  process.  Producing 
stabilized  augmented  forms  is  more  of  an  art  than  a  science.  Although  it  is  true  that  a 
formal  the  process  has  a  formal  analog  in  control  theory,  few  people  working  in  partitioned 
analysis  procedures  have  technical  training  that  would  allow  them  to  exploit  the  analogy. 
The  augmentation  process  is  illustrated  with  examples  in  Section  4.1. 

STEP  6.  Accuracy  Analysis. 

For  a  small  but  important  class  of  coupled  systems  it  is  possible  to  perform  a  fairly  general 
accuracy  analysis  through  the  method  of  the  limit  differential  equation.  For  more  general 
cases  a  Fourier  analysis  technique  can  be  used,  again  with  the  help  of  the  computer  to 
generate  numerical  damping  and  phase  distortion  plots  as  a  function  of  the  temporal 
sampling  rate  [20]. 

Less  thorough  but  more  common  in  practice  is  the  use  of  test  problems  whose  exact 
solution  is  known.  This  has  the  added  advantages  of  verifying  the  implementation  and 
correlating  time  discretization  errors  with  other  modelling  errors. 


SECTION  4 
EXAMPLES 


4.1  STRUCTURE  SUBMERGED  IN  DAAi  FLUID. 


Problem  Description 


A  three-dimensional  structure  is  submerged  in  an  infinite  acoustic  fluid.  A  compressive 
shock  wave  propagates  through  the  fluid  and  impinges  on  the  structure.  The  transient 
response  of  the  latter  is  of  primary  interest. 

The  submerged  structure  is  discretized  by  finite  element  methods.  The  response  of  the 
external  fluid  is  modelled  by  the  first-order  Doubly  Asymptotic  Approximation  (DAAi) 
of  Geers  [10,11].  Through  it  the  fluid  is  replaced  by  a  “membrane”  that  surrounds  the 
structure,  and  which  is  discretized  by  boundary  element  methods. 

The  nature  of  the  interaction  between  the  structure  and  the  DAAi  boundary  may  be 
depicted  as  follows. 
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More  precisely:  at  the  structure  boundary,  the  normal  velocity  of  the  structure  must  match 
the  normal  fluid-particle  velocity  of  the  fluid.*  In  turn,  the  fluid  pressure,  equal  to  incident 
+  scattered,  acts  on  the  structure  as  an  external  force  field. 


Governing  Equations 

For  this  problem  the  global  solution  vector  is  chosen  as 


where  u  is  the  structure  nodal  displacement  vector  and  q  is  the  time-integral  of  scattered 
pressures  ps  («.e.  ps  =  q)  at  the  boundary  element  control  points.  The  governing 
differential  equations  are 
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There  are  no  tangential  velocity  constraints  because  an  acoustic  fluid  is  inviscid. 


where  M„  C4  and  K,  are  the  structural  mass,  damping  and  stiffness  matrices,  respectively, 
N  is  a  pseudo-force  vector  of  structural  nonlinearities,  pj  and  V/  are  vectors  of  incident 
pressure  and  incident  fluid-particle  velocity,  respectively,  at  the  boundary-element  control 
points,  M/  is  the  added  mass  matrix,  A/  is  a  diagonal  matrix  of  boundary-element  areas, 
p  and  c  are  the  fluid  density  and  speed  of  sound,  respectively,  T  is  a  transformation  matrix 
that  related  structural  nodal  forces  to  fluid  forces  at  boundary  element  control  points,  and 
a  prime  denotes  matrix  transposition. 

The  staggered  partition  of  (35)  is 
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there  being  no  need  to  partition  K,  which  is  block  diagonal. 

The  Model  Problem 

Following  a  standard  generalized-coordinate  reduction  process,  the  homogeneous  portion  of 
(35),  with  structural  damping  neglected,  can  be  transformed  to  the  following  dimensionless 
model  problem 


£x  +  w2x  =  — y 
y  +  py  =  x 

Here  x  and  y  represent  generalized  structure  displacement  and  generalized  scattered  fluid 
pressure,  respectively;  £,  w  and  p  are  dimensionless  non-negative  model  parameters.  In 
(36),  dots  denote  derivatives  with  respect  to  a  dimensionless  time  ct/l,  l  being  a  charac¬ 
teristic  problem  dimension  (for  example,  the  radius  of  a  spherical  or  cylindrical  shell) . 

A  thorough  stability  study  of  system  (38)  treated  by  an  implicit-implicit  staggered 
solution  procedure  was  performed;  details  are  given  in  [17].  Only  conditional  stability  was 
attained.  The  largest  stable  stepsize  turned  out  to  be  2 c/((l).  This  was  adjudged  to  be 
unacceptable  small  for  envisioned  applications.  In  fact,  the  limit  is  of  the  order  of  the 
maximum  stepsize  of  an  easily  implementable  fully  explicit  scheme.  Hence,  attention  was 
directed  to  modifying  the  governing  equations  by  augmentation  techniques. 

Augmentation 

An  augmented  system  may  be  easily  obtained  by  proceeding  as  follows.  Solve  the  second 
of  (38)  for  y  =  x  —  py,  replace  this  in  the  right-hand  side  of  the  fir3t  equation,  and  transfer 
the  x  term  to  the  left-hand  side: 

£x  +  x  +  u2x  =  py 
y  +  py  =  x 

The  result  of  the  augmentation  process  is  a  modified  structure  equation,  which  has  acquired 
a  damping  term  x.  This  is  a  radiation  damping  term  resulting  from  the  transportation  of 
high-frequency  kinetic  energy  by  outgoing  waves.  In  the  corresponding  staggered  solution 
procedure,  term  y,  which  is  a  pressure  integral,  has  to  be  predicted;  hence  the  name 
pressure-integral  extrapolation  procedure. 


(39) 


Another  augmented  system  can  be  derived  by  differentiating  the  second  of  (38),  mul¬ 
tiplying  through  by  £,  replacing  £x  =  —  y  —  u2x  obtained  from  the  first  equation,  and 
transferring  the  y  term  to  the  left  hand  side: 

(i  +  z  +  w2z  =  -y 
£y+  (1  +  £fj.y)  =  u2x 

This  is  an  augmented  system  with  a  modified  fluid  equation.  In  the  corresponding  staggered 
solution  procedure  term  x,  which  is  a  structural  displacement,  has  to  be  predicted;  hence 
the  name  displacement  extrapolation  procedure.  More  complex  augmented  forms  may  be 
obtained  by  modifying  both  equations. 

It  is  shown  in  [17]  that  both  augmented  forms  display  significantly  better  stability 
properties,  than  (38).  Unconditional  stability  was  in  fact  attainable  by  judiciously  com¬ 
bining  time  integration  and  predictor  formulas. 

Selecting  an  Augmented  Form 

An  accuracy  analysis  showed  that  the  pressure-integral  extrapolation  procedure  based  on 
(39)  had  generally  better  accuracy  than  the  displacement-extrapolation  procedure  based 
on  (40)  for  the  same  timestep.  However,  the  latter  was  selected  as  the  basis  for  imple¬ 
mentation  because  it  best  satisfied  software  modularity  constraints,  which  called  for  no 
internal  modifications  being  made  to  the  software  component  that  treats  the  structure. 

The  matrix  counterpart  of  (40)  now  forms  the  basis  of  the  USA  (Underwater  Shock 
Analysis)  code  [4].  In  this  code  the  structure  is  integrated  by  the  trapezoidal  rule,  the 
fluid  by  the  3-step  Park  method,  and  a  three-step  predictor  derived  from  a  least-squares 
fit  is  used. 

4.2  STRUCTURE  SUBMERGED  IN  A  C AVITATIN G  FLUID. 

Problem  Description 

A  structure  is  submerged  in  a  fluid  idealized  as  an  infinite  acoustic  medium  incapable 
of  transmitting  tensile  stresses.  A  compressive  shock  wave  propagates  through  the  fluid 
and  impinges  on  the  structure.  If  the  structure  is  sufficiently  flexible  and  the  ambient 
hydrostatic  pressure  sufficiently  low,  the  scattered  negative  pressure  wave  may  induce 
cavitation  in  the  subregion  that  was  traversed  by  the  incident  shock  wave  before  reaching 
the  structure.  This  phenomenon  is  known  as  hull  cavitation. 

Because  of  the  nonlinear  nature  of  cavitation,  a  boundary-element  treatment  of  the 
entire  fluid  domain  as  a  “DAA  membrane”  surrounding  the  structure  is  ruled  out.  (Bound¬ 
ary  element  methods  are  restricted  to  homogeneous  linear  domains.)  Instead,  a  realistic 
computer  analysis  of  this  problem  requires  the  consideration  of  the  three  interacting  fields: 
structure,  cavitating  fluid  volume  and  DAA  membrane.  The  latter  should  be  placed  as  far 
away  as  necessary  to  encompass  the  cavitating  fluid  subregion.* 

*  Inasmuch  as  the  extent  of  the  extent  of  the  cavitation  area  is  generally  unknown  before¬ 
hand,  some  iterations  on  the  placement  of  the  DAA  membrane  may  be  necessary  in  complex 
problems. 


(41) 


The  nature  of  the  interaction  among  the  three  fields  can  be  pictorially  illustrated  as 
follows. 
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Governing  Equations 
The  global  solution  vector  is 


where  u,  is  again  the  vector  of  structural  displacements,  0  is  a  vector  of  displacement 
potentials  at  fluid-volume  nodes,*  and  u/  the  vector  of  DAAi  displacements  at  control 
points  of  the  DAA  membrane.** 

The  fluid  pressure  p  is  related  to  the  displacement  potential  a s  p  —  px  =  rp,  where 

Ph  is  the  hydrostatic  pressure. 

The  governing  semi-discrete  equations  are: 
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where  Q  and  H  are  the  capacitance  and  reactance  matrices  for  the  fluid-volume  mesh; 
these  being  formed  through  a  standard  finite  element  treatment.  Other  terms  are  similar 
to  those  described  in  Section  4.1. 


Stability  Considerations 


Computational  considerations  dictated  that  the  fluid-volume  system  be  processed  explicitly 
using  a  central  difference  formulation.  Hence  conditional  stability  is  accepted  at  the  outset; 


*  Rationale  for  selecting  the  displacement  potential,  which  is  a  scalar  field,  over  a  displacement 
vector  field,  are  elaborated  upon  in  [9]. 


**  The  DAAi  equations  have  been  “turned  around”  to  establish  conjugacy  with  respect  to  the 
fluid-volume  grid. 


%  Vs  s  s  %;.v- 

J  1  J  A. 


LTUAiM 


the  maximum  stepsize  cannot  exceed  the  Courant  limit,  which  is  the  travel  time  of  a  signal 
(shock  wave)  over  the  smallest  dimension  in  the  fluid  mesh. 

The  staggered  solution  method  was  thus  designed  with  the  goal  of  minimizing  the 
reduction  of  the  Courant  limit  due  to  the  presence  of  the  other  fields.  A  stability  analysis 
showed  that  the  presence  of  the  numerical  damping  coefficient  0  could  be  taken  advantage 
of  to  achive  that  goal. 

Accuracy 

The  method  was  implemented  and  stability  predictions  were  verified  on  simple  one¬ 
dimensional  problems.  But  an  accuracy  problem  showed  up:  outgoing  waves  were  not 
properly  radiated  at  the  DAA  boundary,  resulting  in  undesirable  buildup  of  pressure  oscil¬ 
lations.  The  problem  was  solved  by  a  partial  augmentation  of  the  fluid-volume  equations 
with  pc  terms  at  the  DAA  boundary,  which  resulted  in  a  corrected  pressure  calculation 
scheme  at  such  points.  The  oscillations  disapeared. 

Another  type  of  spurious  oscillations  occurs  within  and  near  cavitating  areas;  this 
phenomenon  was  observed  by  Newton  [15],  who  termed  it  frothing.  This  is  controlled 
through  the  numerical  dissipation  term. 

4.3  FLOW  IN  SATURATED  POROUS  MEDIA. 

The  two  fields  are  a  porous  soil  treated  as  an  elastic  media  and  a  percolating  fluid.  The 
interaction  diagram  is  similar  to  that  of  Section  4.1: 


displacements 


Porous 
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Pore 

Fluid 


pressures 


the  main  differences  being  that  displacements  rather  than  velocities  axe  involved,  and  that 
the  interaction  occurs  through  the  volume.  The  latter  observation  influences  computational 
considerations. 

Governing  Equations 

The  semi-discrete  equations  that  govern  the  title  problem  are 

(a  :h:hj  jk:m:  §](;)-{;:)  « 

where  u  and  p  axe  structural  displacements,  p  pore  pressures,  M  and  K  the  soil  mass 
and  stiffness  matrices,  G  and  H  the  pore  capacitance  and  resistance  matrices,  and  M  the 
permeable  mass  matrix. 

A  staggered  solution  procedure  is  proposed,  with  extrapolation  on  the  pore  pressure 
p.  A  stability  analysis  of  this  scheme  show  no  stability  if  M  is  nonzero,  and  conditional 
stability  if  this  term  is  set  to  zero. 
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Augmentation 

This  system  can  be  augmented  in  a  manner  similar  to  that  used  in  the  structure-DAAi 
analysis.  Solve  for  u  from  the  structural  equation,  insert  into  the  differentiated  second 
equation,  and  move  pressure  terms  to  the  left  side.  The  result  is 
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A  stability  analysis  now  reveals  that  the  augmented  equation  treated  by  the  trapezoidal 
rule  and  the  staggered  partition  is  unconditionally  stable  if  M  --  0,  which  is  a  widely  used 
approximation.  If  this  matrix  does  not  vanish,  further  manipulations  are  needed. 

4.4  SUMMARY. 

The  foregoing  examples  have  been  chosen  because  they  represent  working  implementations 
of  partitioned  analysis  procedures. 

The  power  of  the  augmentation  technique  for  stabilization  is  evident  from  the  first 
and  third  examples,  in  which  unconditional  stability  was  the  goal. 

In  the  second  example,  conditional  stability  was  accepted  at  the  outset  because  of  the 
explicit  treatment  of  one  of  the  three  fields.  In  this  case  a  “boundary  augmentation”  was 
beneficial  for  accuracy  reasons. 
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